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Abstract 

The stability of idealized shear flow at long wavelengths is studied in detail. 
A hydrodynamic analysis at the level of the Navier-Stokes equation for small 
shear rates is given to identify the origin and universality of an instability 
at any finite shear rate for sufficiently long wavelength perturbations. The 
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analysis is extended to larger shear rates using a low density model kinetic 
equation. Direct Monte Carlo simulation of this equation is compared with a 
hydrodynamic description including non Newtonian rheological effects. The 
hydrodynamic description of the instability is in good agreement with the 
direct Monte Carlo simulation for t < 50to, where to is the mean free time. 
Longer time simulations up to 2000io are used to identify the asymptotic state 
as a spatially non-uniform quasi-stationary state. Finally, preliminary results 
from molecular dynamics simulation showing the instability are presented and 
discussed. 

PACS numbers: 47.20.Pt, 47.15.Fe, 05.20.Dd, 05.60,+w 
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I. INTRODUCTION 



Uniform shear flow is a prototype non-equilibrium state admitting detailed study at 
both the macroscopic and microscopic levels via theory and computer simulation. This is 
an idealized version of shear flow between parallel plates in which the velocity profile is 
exactly linear in a coordinate orthogonal to the flow direction (as in Couette flow) and 
with a spatially uniform temperature and pressure (in contrast to Couette flow). There 
is a single scalar control parameter, the shear rate a, which measures how far the system 
is driven from equilibrium. This flow is generated by periodic boundary conditions in the 
local Lagrangian frame (Lees-Edwards boundary conditions) that can be implemented at 
the levels of hydrodynamics, kinetic theory, and Newtonian mechanics []TJ— ^] . Although 
these boundary conditions are non-local and therefore not reproducible in real experiments, 
they are ideally suited for computer simulation of this special nonequilibrium state and 
for more penetrating theoretical analysis. In this way, a quantitative study of rheological 
properties usually associated with complex molecular systems has been performed for simple 
atomic fluids 0]. The most complete studies have been via molecular dynamics simulation of 
Newtonian dynamics at high densities and, more recently, by Monte Carlo simulation of the 
Boltzmann equation at low densities HQ. Molecular dynamics simulations have revealed a 
transition from fluid symmetry to an ordered state at sufficiently high shear rates ]7|], which 
has been attributed to a short wavelength hydrodynamic instability ||. The objective here 
is to show that uniform shear flow also is unstable at sufficiently long wavelengths, for any 
finite value of the shear rate [PJTOfl. This instability has not been seen in earlier computer 
simulations due to the small system sizes considered, with consequent restrictions to shorter 
wavelengths. The instability is identified theoretically from a hydrodynamic analysis both 
near and far from equilibrium. This analysis is confirmed quantitatively at short times by 
Monte Carlo simulations of an associated low density kinetic equation. The asymptotic 
evolution of this instability also is explored via Monte Carlo simulation showing transition 
to a non-steady, spatially inhomogeneous state superimposed upon the uniform shear flow. 
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Finally, the instability is demonstrated by new molecular dynamics simulations on larger 
systems to allow consideration of long wavelengths. This instability does not invalidate or 
compromise in any way previous studies of uniform shear flow, since we also establish that 
the flow is stable at sufficiently short wavelengths. 

The boundary conditions generate viscous heating so that uniform shear flow is not 
stationary. This viscous heating can be controlled by the introduction of a non- conservative 
external force that acts as a uniform thermostat. The simplest choice is a Stokes law drag 
force on each particle proportional to its velocity relative to the local macroscopic flow. 
The proportionality "constant" then is adjusted to compensate for the heating. There are 
several possibilities in the detailed implementation of the thermostat, leading to the same 
properties for the stationary state but different hydrodynamics for small perturbations from 
that state. The theory and Monte Carlo simulations are carried out for both global and 
local thermostats. The role of the thermostat is studied in the Appendix where it is shown 
that the qualitative features of the instability are not sensitive to the choice of thermostat. 

In the next section, the usual Navier-Stokes hydrodynamic equations are considered. 
These equations are restricted to small gradients of the hydrodynamic fields relative to 
equilibrium, and consequently the shear rate must be small in this analysis. Otherwise, 
the density and interatomic force law can be considered arbitrary within the fluid phase. 
The stationary solution for uniform shear flow is identified, and the linear hydrodynamic 
equations for small perturbations of this solution are studied. The five hydrodynamic modes 
are determined in detail for the special case of spatial perturbations orthogonal to the 
flow. A critical wavevector, k c (a), is determined such that for wavevectors k < k c (a) the 
perturbations grow as a function of time. The critical wavevector vanishes as the shear 
rate a goes to zero, but for any finite value of the shear rate there are sufficiently small 
wavevectors (long wavelengths) such that the perturbation is unstable. 

These hydrodynamic predictions are tested by comparison to Monte Carlo simulations at 
the more fundamental kinetic theory level. A model kinetic equation has been analyzed for 



states near uniform shear flow, without restriction on the shear rate |I| . The hydrodynamic 



equations for small deviations from uniform shear flow determine the critical wavevector, 
k c (a), for values of the shear rate a beyond the limitations of the Navier-Stokes equations 
where efficient Monte Carlo simulations are possible. The theoretical prediction of the 
growth of initial perturbations is compared with a direct Monte Carlo simulation of a 
solution to the kinetic equation. The results confirm both the hydrodynamic analysis and 
the prediction of an instability for times up to about 50to, where to is the mean free time, 
after which the initial growth has exceeded the limitations of the linear stability analysis. 
The simulation results are continued up to 2000to to explore the ultimate stabilization by 
non-linear effects. The asymptotic state for the hydrodynamic fields appears to be a system 
size dependent standing wave with a period of about 50t - Further details and discussion 
are given in section |T|. 

In section |Iy| some preliminary attempts to see the instability at high densities using 
molecular dynamics for the hard sphere fluid are discussed. The system dimension in the 
direction of the spatial perturbation is increased by an order of magnitude relative to previous 
simulations. A long wavelength perturbation is found to grow on a very long simulation time, 
with no indication of approach to the steady uniform shear flow. A quantitative comparison 
with theory at the required larger densities and shear rates is now possible using a recently 
developed kinetic model for the hard sphere Enskog equation for application although 
the details have not been worked out at this point. The results of the theory and simulations 
are summarized and discussed in section [V|. 



II. NAVIER-STOKES ANALYSIS 

On sufficiently large space and time scales the dynamics of a fluid is well-described by 
hydrodynamic equations obtained from the exact conservation laws for the average mass, 
energy, and momentum densities, together with approximate constitutive equations for the 
associated fluxes. For the analysis here we use the density, n(r,t), temperature, T(r, t), and 
flow velocity, U(r, t), as dependent variables rather than the density, energy density, and 
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momentum. The general form of these conservation laws is 

D t n + nV • U = 0, 







(2.1) 



D t T + (7 - l)a _i V ■ U + (mnC v )- L (V ■ q + P i] d ] U l -w) = 0, (2.2) 

D t Ui + p- l d iP + p-'^Py = 0, (2.3) 

where D t = dt +U- V is the material derivative. The parameters occurring in these equations 
are the mass density, p = mn, the specific heat at constant volume, C v , the pressure, p, the 
ratio of specific heats at constant pressure and volume, 7 = C p /C v , and the coefficient 
of expansion, a = —rT l (dn/dT) . These parameters are the same functions of the local 
density and temperature as in equilibrium. Finally, the irreversible heat and momentum 
fluxes are denoted by q and Pij, respectively, and w = w(n,T) is the rate at which work is 
done by the external force representing the thermostat. Its detailed form will not be required 
in this section. 

The above equations are incomplete until constitutive equations for the fluxes are spec- 
ified in terms of the hydrodynamic fields. However, the special solution of uniform shear 
flow exists independent of this choice. It is defined by a spatially constant temperature, 
density, heat flux and momentum flux, and a flow velocity whose only non-vanishing com- 
ponent is U ajX = ay. The shear rate, a, provides the single control parameter measuring the 
deviation from equilibrium. The boundary conditions are simple periodic conditions in the 
local Lagrangian coordinate frame, r' = r — U s (r)i. Substitution of these assumptions for 



n = n s , T = T s , and U = U s into the above conservation laws shows that (|2.1|) and ( [2.3| ) 
are satisfied, while ( |2.2| ) reduces to 

d t T s = -(mrisC^s)' 1 (aP S)Xy - w(n s , T s )) . (2.4) 

This expresses the temperature evolution as a competition between the viscous heating 
effect, oc aP S)Xy , and the cooling by the thermostat, oc w(n s , T s ). A steady state is obtained 
by choosing the thermostat to cancel the viscous heating, 
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aPs,xy = w(n s ,T s ). (2.5) 

The various thermostats described in the Appendix all satisfy ( |2.5| ) but differ for states away 
from the steady state. 

Next consider the equations for small deviations of the hydrodynamic variables from 
the uniform shear flow state, retaining only terms linear in these deviations. To proceed 
it is necessary to specify the constitutive equations for the heat and momentum fluxes. In 
this section, attention is limited to small spatial gradients, including the shear rate, so that 
Fourier's law and Newton's viscosity law apply, 

q = -kVT, (2.6) 

Pij = -v(d i U J + djU, - ■ U) - rfSijV ■ U. (2.7) 

Here re(n, T) is the thermal conductivity, 77(71, T) is the shear viscosity, and rj'(n, T) is the 
bulk viscosity. It follows immediately that 

q s = 0, P s>ij = -rj s a (5 ix 5 jy + 5i y 5 jx ) , (2.8) 
<5q = -k s V5T, (2.9) 

2 

8 Pij = -rj s (di5Uj + dj5Ui - -%V ■ 5\J) - r/^%V ■ 5XJ-a (5 ix 5 jy + 5 iy 5 jx ) (ri s>n 5n + r] S)T 5T) . 

(2.10) 

An abbreviated notation has been introduced where the subscript s on a quantity in- 
dicates it is evaluated at n s ,T s . Also, p S)Tl = dp(n s ,T s )/dn s , p St T = dp(n s ,T s )/dT s , 
7/ Sjn = dr](n s ,T s )/dn s , t] s ^t = dr](n s ,T s )/dT s , etc. With these results, the closed set of 
linear hydrodynamic equations for perturbations of uniform shear flow at small shear rates 
are given by 

(d t + V s - V)6n + n s V -6V = 0, (2.11) 
7 



(d t + U s • V) 5T + ( 7s - l^V • 5V + K^J- 1 (-k s V 2 5T - 2 Vs a (dJU y + d y 5U x )) 

= (mrisC^s)" 1 (a 2 {i] S)n 5n + r] SjT ST) + 5w^ , (2.12) 

(dt + U s • V) SUi + 5 ix a5U y + pj 1 (ps^diSn + Ps , T di5T) + p^d.SP^ = 0. (2.13) 

In this section we choose a local thermostat for which 5w in the temperature equation 
compensates for the excess viscous heating due to perturbations of the temperature and 
density, 

5w = -a 2 (r] s ,Jn + r] S:T 5T) . (2. 14) 

This does not imply a constant local temperature, however, except for spatially homogeneous 
deviations from uniform shear flow, although it does lead to a constant average temperature 
or kinetic energy for the whole system. 

These differential equations have constant coefficients which suggests an equivalent al- 
gebraic form using Fourier and Laplace transformation. Consider first the Fourier trans- 
form. The boundary conditions are periodic in the local Lagrangian frame given by 
r[ = Ti — U S j(r)t = Aij(t)rj, where Ajj(t) = 5ij — aSi X Sj y t, so it is appropriate to define 
the transform with respect to the variable r', 

Sy a (k,t) = J dv'^ r '5y a {v,t) = J dr^- r 5y a ( r ,t), (2.15) 

where 5y a (r,t) denotes the set of perturbations, considered as a function r' in the first 
equality. The periodic boundary conditions require fc, = 2nj7r/Lj, where rii are integers 
and Li are the linear dimensions of the system. However, since the time derivative in the 
hydrodynamic equations is taken at constant r, the representation following the second 
equality is useful with ki(t) = kjAj^t). The Fourier transformed hydrodynamic equations 
are 

d t 5y a + (A(a) - ik(t)B(a) + k 2 (t)D)^ 5y p = 0. (2.16) 
The three matrices A (a), B(a), and D are 
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(2.17) 



B a p{a) 













n s k x n s ky Ti s k z 

c\k x — 2ai] s C2k y c\k y — 2ar) s cik x c\k z 



\ 



Ps 1 (k x Ps,n ~ k y aris,n) Ps 1 (k x p s ,T - k y ar] SjT ) 
Ps 1 (kyP s , n - k x ar] s ^ pj 1 (k y p s , T - k x ar] SjT ) 

Ps^zPs^ P^zPs,T 





















(2.18) 



Dafl = P. 



-1 








c s h x h z 



(2.19) 





p s c 2 K, s 

a s k\ + i] s a s k x k y 

o s k y k x a s k 2 y + i] s a s k y k z 

^ <r s k z k x cr s k z ky Gsk\ + f] s ^ 

where k = k(t) is the unit vector along k(t), c\ = (7^ — l)a7\ c 2 
a s = lrj s + r)' s . 

To simplify the analysis attention is restricted here and below to spatial perturbations 
only along the velocity gradient direction, i.e. k =k y. In this case the linear hydrodynamic 
equations have time-independent coefficients (i.e., k(t) = k), 



(mn s C VyS ) \ and 



d t 5y a + F a/3 (k, a)5y p = 0, F af} {k, a) = (A{a) - ikB(a) + k 2 D) 



a/3 



(2.20) 



and the matrices B(a) and D simplify to 



B a/3 (a) 














n s 













-2af] s c 2 


Cl 







-Ps 















Ps 































(2.21) 
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C 2 K S 

o o p -% 
o o 

















Ps'ilVs + v's) o 
v o pj% j 

Equation (|2.20|) can be solved by Laplace transformation, 



(2.22) 



[k,z) = / dte- tz 5y a (k,t), 
Jo 



(2.23) 



with the result 



y a (k, z) = [zl + F(k, a)U5y p (k, t = 0). 



(2.24) 



The eigenvalues oj^(k, a) of the matrix F(k, a) define the five simple hydrodynamic poles at 
z = -cu w (k,a). The resulting five exponentials in time represent the hydrodynamic modes 
for relaxation of the perturbations around uniform shear flow. At a = (perturbations of 
equilibrium) they are the two sound modes, a heat mode, and a two-fold degenerate shear 
mode. For finite shear rate, the modes are more complicated and have qualitative differences. 
To illustrate, consider first the case of k — > at fixed, finite a, 



CO 



W(k,o) 



hk 2 



\ 



-|(1 + i V3)b 2 (a)k 2 / 3 + |(1 - i v / 3)6 3 (a)A; 4/3 + b A k 2 
-i(l - i y/3)b 2 (a)k 2 / 3 + + i V3)b 3 (a)k^ 3 + b 4 k 2 
b 2 (a)k 2 / 3 + 6 3 (a) A; 4 / 3 + b 4 k 2 
(Vs/ps) k 2 



(2.25) 



The coefficients in these expressions are real, 



b\ = (Vs,TPs,n ~ Vs,nPs,T) /mp s ,T, b 2 (a) = (2a 2 r] s p s>T / p 2 s C Vy 



1/3 



63(a) = (2a 2 7] s ri^ T {p s C VtS ) 1 + n s p s ^ n + (j s - 1) a s x p„t) A° s b 2 (a), 
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h = i (- h + k s (psC^y 1 + p- 1 (2 Vs + a s )) . (2.26) 

There are two diffusive modes, ~ k 2 , but the other modes are non-analytic about k = and 
represent more complex spatial dependence. This behavior can be traced to the fact that 
the matrix A(a) — ikB(a) is not normal and cannot be diagonalized. Thus, at fixed a ^ 
there is a crossover in the transformation of F a ^(k, a) from normal diagonal form to Jordan 
form at small k. This is reflected in the eigenvalues if they are expanded in k at fixed a, 
as above. For similar reasons, care must be used in representing Sy a (k, t) as an expansion 
in the hydrodynamic modes with constant coefficients, since at small wavevectors there is a 
crossover to a mode expansion whose coefficients have algebraic time dependence. 

The two propagating modes in Q2.25 ) are unstable, since 62(a) > 0. The above Navier- 



Stokes analysis applies for small but finite shear rate, and small k (long wavelength). The 
expansion in k verifies that the asymptotic long wavelength modes are always unstable. 
At larger values of k the modes are again stable, as follows from an exact evaluation of 
the eigenvalues. There is a critical wavevector, k c (a), such that for k > k c (a) the modes 
are stable whereas they are unstable otherwise. These qualitative results apply without 
restriction to the atomic force law, density, or temperature. Figure [I] shows k c (a) as a 
function of a for the special case of hard spheres at three densities, n* = na 3 = 0.0, 0.2, 
and 0.4 (k and a measured in units of the inverse mean free path and mean free time for 
the hard sphere Boltzmann equation). The thermodynamic properties are calculated using 
the Percus-Yevick approximation, while the transport coefficients are calculated from the 
Enskog kinetic theory. 

The instability is due to three matrix elements, -B42, -B23, an d and is therefore 
present at order k. The density perturbation is constant to this order and we choose Sn = 
to simplify the discussion. The relevant variables are then the temperature perturbation 
ST, the longitudinal velocity perturbation SU y , and the transverse velocity perturbation SU X 
which to this order obey the equations 

8 t ST + ( 7s - l)a- l d y SU y - 2a Vs (psC^)' 1 d y SU x = 0, (2.27) 
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d t 8U x + a5U y - p- x afis iT d v 5T = 0, d t 5U y + p^p^dyST = 0. (2.28) 

The first equation has a coupling to the transverse velocity field due to the reference shear 
flow; the second equation provides a feedback to the temperature equation through the 
same shear flow. These couplings alone would lead simply to a renormalization of the 
sound velocity. However, the shear flow also couples the transverse field to the longitudinal 
field for an additional feedback mechanism to the temperature equation through the pressure 
gradient. This second mechanism is responsible for the instability. For very long wavelengths 
the above equations can be simplified to give 

dfST ~ (2a 2 r lsPs , T /p 2 s C v , s ) d 2 5T, (2.29) 

which exhibits the unstable modes of ( |2.25| ). 



III. MONTE CARLO SIMULATION 

The hydrodynamic description can be derived from a more fundamental level of kinetic 
theory. In principle, this also allows derivation of hydrodynamic equations without the 
restriction of the Navier-Stokes approximation to small shear rates. A model kinetic theory 
for the practical calculation of such generalized hydrodynamic equations is given in Ref. 
[fToH . The resulting equations are limited to long wavelengths, as in the Navier-Stokes case, 
but the reference state of uniform shear flow can have a large shear rate. An analysis 
of the hydrodynamic modes shows that there is a critical wavevector similar to that of 
Figure |l|. It is possible to test this hydrodynamic description by a direct simulation of the 
more fundamental solution to the kinetic equation. For practical reasons the simulation is 
more efficient at larger wavevectors and shear rates than can be justified by Navier-Stokes 
hydrodynamics, and this is the primary reason for considering the more complex generalized 
hydrodynamics. 

The kinetic equation is a single relaxation time BGK equation [13 given by 

(Jf t + v • V r + V.m- 1 • F cx ^j /(r, v, t) = -v(f(r, v, t) - £(r, v, *)). (3.1) 
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Here F ext is the external force representing the thermostat, fe is the local equilibrium dis- 
tribution, and v is an average collision rate. The exact stationary solution for uniform 



shear flow in the presence of a thermostat has been studied in detail [0. A variant of 
the Chapman-Enskog method can be used to study normal solutions for states near uniform 
shear flow, and to obtain the corresponding hydrodynamic equations for small perturbations 
relative to this state. The complete details of this solution and the hydrodynamic modes as 
a function of the shear rate can be found in [ lOfl . The method of simulation is a variant of the 



Bird Direct Simulation method [15] . This is a Monte Carlo technique composed of two parts 



at each discrete time step: free streaming and collision. The volume of the system is divided 
into cells of dimension smaller than the mean free path, and iV particles are introduced at 
t = with positions and velocities sampled statistically from a specified initial distribution 
function. The distribution of particles is calculated at t = At, with At much smaller than 
the mean free time, as follows. First the particles are displaced by {Ar M = v M At}. Next 
the velocity of each particle fi is replaced with probability z/(n M , T M ) At by a random velocity 
sampled from the local equilibrium distribution, fe({n^, T M , U M }; v). Here n^, T M , and U M 
are the density, temperature, and flow velocity in the cell containing particle /i. In this col- 
lision stage, strict conservation of momentum and energy may be violated due to statistical 
fluctuations. To compensate for this artificial effect, the velocities of the particles in each 
cell are conveniently displaced and rescaled. The whole process is then repeated for each 
subsequent time step. In this way, the one particle distribution function /(r,v,t) (coarse 
grained over the cells) is determined, from which the hydrodynamic fields can be computed 
directly as averages. 

In our simulations we have considered a system of size L = 2n/k, with k = O.l(foto) 1 ! 
along the ^/-direction at a shear rate a = OMq 1 . Here to = l/K^s) is the mean free 
time with v oc n (Maxwell molecules) and v = (2kBT s /m) 1 ^ 2 is the thermal velocity. In 
the remainder of this section, we take to — 1, Vq — 1, T s — 1, and n s — 1. Lees-Edwards 
boundary conditions are used to drive the shear flow [|T|. These are simple periodic boundary 
conditions on both the position and velocity variables in the local Lagrangian frame at 
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y = ±L/2. Both local and global thermostats have been studied. The local thermostat 
is the same as that of section |I| except that the pressure tensor is no longer limited to its 
Navier-Stokes form. In addition, two global thermostats are considered for more efficient 
implementation of the simulation. These are described in the Appendix. The analysis there 
and the results of the simulation show that the instability is not sensitive to the choice of 
thermostat. Consequently, only the results using the global thermostats are presented here. 
Starting from a distribution corresponding to uniform shear flow ||, the initial condition 
has been prepared by displacing and rescaling velocities so that 

SU XtV (y, 0) = -SU x , y (0) sin ky, 6U z (y, 0) = 5n(y, 0) = 5T(y, 0) = 0, (3.2) 

with 5U y (0) = 0.1 and SU x (0) = —0.03. The technical parameters of the simulations are 
N = 628 000 particles, a time step At = 0.02, and a cell width AL = 0.05. The data have 
been averaged over 10 different realizations in the simulations of Figs. 0, (|, |^, and [j], and 
over 2 realizations in those of Figs. [5] and |5|. 

First we consider the short time dependence at a fixed position y = — L/4. Figure |2] shows 
5U x (t) = 5U x (—L/4:, t) and 5U y {t) = 5U y (—L/4:, t). The dashed lines are the results from the 
hydrodynamic analysis of the BGK model near uniform shear flow [nj. The good agreement 
up to t ~ 50 shows that the instability is not just a consequence of the assumptions behind 
the hydrodynamic description. The subsequent differences between simulation and theory 
are due to the fact that the latter is limited to small deviations from uniform shear flow. For 
longer times we find both large amplitudes for Sy a and large deviations of the distribution 
function from that of the unperturbed state. 

To investigate the asymptotic state the system, we have performed the simulations for 
much longer times. Figure ^ shows SU x (t) and T(t) = T{— L/4, t) for time up to t = 2000. 
Both the velocity and the temperature oscillate in time and are modulated by a slowly 
varying amplitude relative to their asymptotic average values. The maximum value T w 
9 corresponds to k pa 0.3, which is close to the value k c (a) at which the uniform shear 
flow with a = 0.5 would become marginally stable |§. Thus the initial dynamics tends 
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toward stabilization but does not ever cross over into the stable domain. Consequently, 
the asymptotic dynamics is not simply that of stationary uniform shear flow at a different 
temperature, but rather a quasi-stationary state with different spatial structure. 

The thermostat used in the above simulations allows a global change in the average 
temperature, which is responsible for the amplitude modulations at long times. In order to 
have a more controlled asymptotic state, we have considered a variation of this thermostat 
that maintains the average temperature constant (see the Appendix). The quantities SU y (t) 
and n(t) = n(—L/4,t) are plotted in Fig. |4|. After a transient period of length t ~ 100, a 
stable oscillatory behavior of the velocity appears with a period r ~ 54. The shape of 5y a 
over one cycle is shown in Fig. |5|. Here t' = t — t tr , where t tT = 463.4, and the curves are 
averages over 20 successive cycles. The value of t tT has been chosen as to assure that the 
transient time is over and also with the criterion that SU X = at t' = 0. Inspection of the 
results shows several regularities. First, the following symmetry relation appears: 

5y a (t , + T/2) = ±5y a (t'), (3.3) 

where the minus sign applies to the velocity and the plus sign applies to the density and 
the temperature. Next, at times t' ~ 0.36r, 0.86r, where 5U X have extrema, 8U y , 5T, and 
Sn seem to have nodes. Note also that the vector 5XJ rotates anti-clockwise and that most 
of the time 5T > is correlated to 5n > and 5T < is correlated to 5n < 0. 

The above results indicate the wave character of the asymptotic state. To confirm this 
point, we have analyzed the profiles of the hydrodynamic quantities at relevant times. The 
results are consistent with two independent invariance relations: 

8y a (y,t') = 8y a (y + L/2,t' + r/2), (3.4) 

8va(y, = ±8y a (-y, 0- ( 3 - 5 ) 

Their combination yields 8y a (y,t') = ±5y a (—y — L/2,t f + r/2), which implies Eq. 
( 3.3[) . Figure | shows the spatial variation of 8n(y,t') and 8T(y,t') at f = 

15 



0, 0.14r, 0.25r, 0.36r, and 0.5r. Not shown are times 0.5r < t' < t because they can be 
reproduced by use of the relation ( |3.4j ). As observed in Fig. |5| for the special point y = — L/4 
it is seen that a high (low) temperature is generally correlated to a high (low) density. The 
spatial distribution of the temperature is highly non-uniform even though the thermostat 
maintains a constant average temperature. Figure |7] shows a vector representation of the 
components of 5XJ throughout the system at the values of t' . As anticipated from Fig. || 5XJ 
rotates anti-clockwise throughout the system. The layers y — 0, ±L/2 are always nodes of 
5XJ and extrema of 6T and 5n. The pattern indicated by Figs. |6| and [7| can be described as a 
periodic standing wave represented by the superposition of two symmetrical waves travelling 
in opposite directions. 

In summary, for initial values of k and a in the predicted unstable domain small pertur- 
bations of the hydrodynamic fields grow according to the linear hydrodynamic equations for 
t < 50. Subsequently, non-linear effects invalidate this theoretical analysis. The simulations 
show a transient period up to about t « 100, after which a quasi-stationary state is observed 
for 100 < t < 2000. In this asymptotic state the vector quantities oscillate at a period ap- 
proximately twice that of the scalar fields. The oscillations are spatially non-uniform for all 
fields considered. 



IV. MOLECULAR DYNAMICS SIMULATION 

The most extensive prior studies of uniform shear flow have been for dense systems 
via molecular dynamics (MD) simulations. No signature of the instability discussed here 
has been noted in these previous results. There are two significant differences between 
the theory and simulation discussed above and the detailed implementation of earlier MD 
simulations. The first, and perhaps most important, has been the consideration of small 
system sizes relative to the wavelengths necessary to see the instability. Typically, system 
size is determined by the simulation time such that a sound wave will not traverse the 
system and generate correlations. At high densities this has led to consideration of system 
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sizes small compared to the critical wavevector for instability. A second difference from the 
discussion above is the method for imposing a thermostat. In MD simulations it is efficient 
to impose the temperature control by a global re-scaling of the velocities only after as many 
as 100 collisions. However, at the shear rates considered this implies significant heating 
between applications of the thermstat and the temperature is more like a sawtooth in time 
rather than constant. Thus it is difficult to give a direct theoretical correspondence to the 
MD simulation results. 

To demonstrate the existence of the instability using MD simulations we have considered 
a rectangular unit cell with one side expanded 15 times larger than the other two ( 1620 hard 
spheres of diameter a at a density of na 3 = 0.5 in a cell of size 6a x 90<j x 6a). The Lees- 
Edwards boundary conditions are applied on the x, z surfaces so that the velocity gradient is 
along the larger dimension of the system, allowing study of much longer wavelengths in the 
direction of the gradient than have previously been considered. A second difference from 
previous simulations is a more frequent application of the velocity scaling to control the 
themperature and better represent continuous cooling. In our simulations the velocities are 
rescaled whenever the temperature differs by 5% from its set value. The simulation cell is 
divided into a fixed number of subcells, the local velocity field in each subcell is calculated 
and the excess kinetic energy computed relative to the local velocity field. The velocities 
of the particles in each subcell are then rescaled so that the total excess kinetic energy 
of each subcell is equal to the set value. The amount of rescaling is determined locally, 
and thus corresponds to a local thermostat, A (r, t) = A [n (r, t) , T (r, t)], as discussed in 
the Appendix. Our implementation follows Hess [f|] in that we assume uniformity of all 
quantities in the directions perpendicular to the velocity gradient so that the subcells are 
thin slabs and the only spatial variations are in the direction of the velocity gradient, as 
in the Monte Carlo simulations. Due to the large number of particles in the simulations, 
the length of the simulations is relatively modest t « 5 x 10 5 collisions. The shear rate was 
fixed at a = 1 . 77 \JksT /ma 2 and an initial perturbation with k y = 2ir/Ly was monitored. 
Our theoretical estimates indicate this should correspond to conditions of instability. Figure 
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|8] shows the x-component of the velocity field growing steadily throughout the simulation, 
clearly indicating the instability Conversely, perturbations with a wavelength one quarter 
of this value appear to be stable as expected. Similar results are observed for the density 
field as well. Our primary conclusion from these preliminary results is that the instability 
can be observed and studied by MD simulations if larger system sizes are considered and 
more care is taken with application of the thermostat. 

V. DISCUSSION 

Uniform shear flow has been a prototype state for the study of fluids far from equilibrium, 
using both theoretical and computer simulation methods. Until recently, it has been assumed 
that this state is stable except at high densities and short wavelengths. The new results 
reported here and in 0,|1(J show that this simple macroscopic state is unstable at sufficiently 
large wavelengths. Previous studies via simulation have not seen this instability due to 
finite system sizes. However, theoretical analysis at the hydrodynamic level clearly shows 
the mechanism and parameter space for this instability. In the present work this theoretical 
analysis is tested both qualitatively and quantitatively. At the qualitative level, both Monte 
Carlo simulations of a kinetic theory description and molecular dynamics simulation of the 
Newtonian dynamics show clearly that this flow is unstable at long wavelengths. The Monte 
Carlo simulation also confirms quantitatively the predictions of the theory on the time scale 
for which the linear analysis is valid. At longer times the Monte Carlo simulation shows 
clearly that the asymptotic state is spatially non-uniform with a periodic variation in time. 
A precise theoretical description of this final state has not been developed at this time. 

The analysis here has been limited to spatial perturbations along the gradient of the 
velocity in the stationary state. More general perturbations will lead to more complex flows 
due to the coupling to the convective flow. A theoretical analysis of this more general case 
is in progress but is significantly more complex. Due to the fact that the k = matrix of the 
hydrodynamic equations is non-diagonalizable (non-normal), the corresponding eigenvalues 
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do not fix a unique eigenspace. The resulting dynamics is not simply a superposition of 
decaying or growing modes, but rather includes as well algebraic growth. The resulting 
analysis of conditions for instability is more complex and will be reported elsewhere. 
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APPENDIX: ROLE OF THE THERMOSTAT 

The analysis of section |I| made use of a specific choice for the thermostat. Other choices 
are possible and more convenient for computer simulation. In this work we use two different 
types of thermostats, both obtained from an external force at the microscopic level of the 
form 

F cxt (r,t) = -mA(n(r,t),T(r,t))[v-U(r,t)]. (Al) 



The corresponding source term w in the temperature equation (|2.2j ) is 



w{n(r, t),T(r, t)) = -2A>(r, t), T(r, t))A(n(r, t),T(r, t)), (A2) 

where K(n(r, t), T(r, £)) is the kinetic energy density. The thermostat parameter, 
X(n(r, t), T(r, £)), is always chosen to ensure constant temperature in the reference state, 
X s = X(n s , T s ) = —aP SjXy /2K s , as follows from (|2.5| ). In this work we use three different ther- 
mostats, one local and two global. The local thermostat is defined by X(n(r, t), T(r, t)) — > 
A s + X s , n 5n(r, t) + X St TST(r,t), to linear order, so that 

Wl (n(r, t), T(r, *)) -> -2K S X S - 2 (KX) s n 5n(r, t) - 2 (KX) s T 6T(r, t). (A3) 
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The coefficients A s>n and A Sj t are chosen such that all viscous heating proportional to 8n(r, t) 
and 5T(r,t) is compensated by the source term. While there is still a local temperature 
variation, this thermostat holds the average temperature constant. It is this local thermostat 
that is used in section |T], with P xy approximated by its Navier-Stokes form. The other two 
thermostats considered are global, i.e. the parameter A is spatially constant. The simplest 
case is A = A s with a corresponding source term 

w 2 (n(r, t),T(r, t)) -> -2K S \ S - 2\ s K s ,Jn(r, t) - 2X s K s>T 5T(r, t). (A4) 

In this case there are greater effects due to viscous heating by the perturbations than with the 
choice Wi, and even the average temperature of the system changes except at the stationary 
state. A second global thermostat eliminates this global change using the form, A — > X s + 
K,n8n + ^s,t8T, where Sn and 5T are the volume averaged density and temperature. The 
source term in this case is 

w 3 {n{r, t),T(r, t)) -> w 2 (n(r, t),T(r, t)) - 2K S [X s>n 5n + A SjT 5T] . (A5) 

In Fourier representation the sources w 2 and w 3 are the same except for the k = dynamics. 

Both theory and the Monte Carlo simulations of section |TTT| indicate that the details of 
the perturbed dynamics depend on the thermostat used, but that the mechanisim for the 
instability is insensitive to the choice of thermostat. To illustrate this we repeat the analysis 
of section |I| using the global thermostat ( |A4|) (again in the Navier-Stokes limit). There 
are two new terms in the temperature equation. One is due to temperature perturbations 
and leads to viscous heating even for uniform perturbations. The second is due to density 
perturbations and provides a direct k - independent coupling of the temperature and density 
equations (at zero shear rate the density and temperature are only coupled indirectly via 
the gradients of the longitudinal velocity component and the pressure). The matrices B(a) 
and D are unchanged from section [H|, but these new couplings replace A21 and A22 in ( |2.17|) 
with non-zero values proportional to a 2 , 
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A 21 {a) 



\ ( 



-c 2 (a 2 ?7 Sin + w n ) 



(A6) 



y ^22(0) J y-c 2 (a 2 T] StT + w T ) , 

where w n = dw(n s ,T s )/dn s , wt = dw(n s ,T s )/dT s . The matrix A(a) — ikB(a) is now 
diagonalizable and an expansion of the hydrodynamic modes in powers of k rather than k 2 ^ 3 
is obtained, 
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«(k,a) 



A 22 (a) + di(a)k 2 
d 2 (a)k 2 

ic(a)k — (d 3 (a)/a 2 — d<±(a)) k 2 
—ic(a)k — (d 3 (a)/a 2 — d 4 (a)) k 2 
(Vs/Ps) k 2 



\ 



with the coefficients, 



, , s , 2a 2 c 2 ri s 7] sT + c x p SiT p StT (n s r] s ^ n a 2 + n s w n + 2r} s a 2 
dAa) = k s c 2 H 7 ^ ; r~ + 



P s c 2 (?7 S)T a 2 + ty T ) 



P s c 2 (r]s,Ta 2 + w T ) 2 



(A7) 



w» [^sPs.n (a 2 ^,T + u> r ) - P S ,T (a 2 n s rj Sin + n s w n + 2r] s a 2 )] ' 



c 2 (a) 



[w B p B ,n (q 2 ^,t + wt) - Ps,t (a 2 nsVs,n + n s w n + 2r] s a 2 )} 
Ps (a 2 Vs,r + w T ) 



d*(a) 



2p s c 2 (a 2 r? SjT + w T ) 



ci + 



(a 2 n s r] S)n + n s w n + 2r] s a 2 ) 
{a 2 r] S)T + w T ) 



M a ) = ~o rf 2(a) - - ( , 

2 p s (a + w T 



2 Ps 



(A8) 



There are two diffusive modes, two propagating sound modes, and a time-modulated 
diffusive mode. The sound modes are unstable at long wavelengths for shear rates satisfying 
d^{a) — a 2 di{a) > 0. This possibility has been explored in reference || for the special 
case of hard spheres. In that case both d^(a) and d^a) are positive and the modes are 
unstable for sufficiently small shear rates. The corresponding critical wavevector, k c (a), for 
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this thermostat is determined from the exact eigenvalues and shown in Figure |9] for the same 
densities as in Figure |I[ The domain of instability at long wavelengths is similar to that of 
Figure [I] using the thermostat W\. 

The thermostat of this Appendix allows greater viscous heating than that of section [IT] 
for states perturbed from uniform shear flow. The hydrodynamic modes are quite different, 
reflecting a sensitivity to the choice of thermostat. These qualitative differences can be traced 
to the mathematical differences between diagonalizable and non-diagonalizable matrices, 
A a p. Nevertheless, the mechanism for the instability described at the end of section [II 
remains effective in both cases. These conclusions are not limited to the Navier-Stokes 
approximation, but are confirmed as well for the kinetic model results described in section 
I fj for larger shear rates. 
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FIGURES 

FIG. 1. Critical lines for stability for hard spheres at n* = 0.0 (solid curve), n* 
curve), and n* = 0.4 (dashed-dotted curve) with the local thermostat. 



= 0.2 (dashed 



FIG. 2. Plot of 5U x (t) and 5U y (t). The solid lines are simulation results and the dashed lines 
correspond to the analysis of Ref. |l0| . 

FIG. 3. Plot of 5U x (t) and ST (t). 

FIG. 4. Plot of SU y (t) and 8n(t). 



FIG. 5. Plot of 5y a (t') over one cycle. 

FIG. 6. T(y,t') and n(y,t') as a function of y for (from top to bottom) 

t' = 0, 0.14t, 0.25r, 0.36r, and 0.5r. 



FIG. 7. Vector plot representing 5U x {y, t') and SU y (y, t') at t' = 0, 0.14r, 0.25r, 0.36r, and 0.5r. 

FIG. 8. The long-wavelength component of the velocity in the flow direction as a function of 
time from molecular dynamics simulation. Units are m = a = ksT = 1. 

FIG. 9. Critical lines for stability for hard spheres for the same densities as in Fig. || with the 



global thermostat (|A4] ). 
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